Friedel Oscillations疑问

问题

动量空间中电势为

\begin{align} \varphi (\vec{q})= \frac{\varphi_{ext(\vec{q})}}{\varepsilon(\vec{q})} \end{align}

其中

\begin{align} \varphi_{ext}(\vec{q}) = \frac{-e}{\varepsilon_0 V q^2 } \end{align}

所以

\begin{equation} \label{eq:1} \varphi (\vec{q})= \frac{-e}{\varepsilon_0 V q^2 {\varepsilon(\vec{q})}} \end{equation}

Thomas-Fermi 近似的结果

Thomas-Fermi 近似给出

\begin{align} \varepsilon_{TF}(\vec{q}) = 1 +\frac{q_{TF}^2}{q^2} \end{align}

将上式代入 \(\varphi(\vec{q})\) 有

\begin{equation} \varphi (\vec{q})= \frac{-e}{\varepsilon_0 V q^2 \left(1 +\frac{q_{TF}^2}{q^2}\right)} =\frac{-e}{\varepsilon_0 V \left(q^2 +q_{TF}^2\right)} \end{equation}

对其 Fourier Transform

\begin{align} \varphi(\vec{r}) =& \int_{-\infty}^{+\infty} e^{i \vec{q}\cdot \vec{r}} \varphi (\vec{q}) \mathrm{d}^3 q \\ =& \frac{-e}{2 \pi^2 \varepsilon_0} \cdot \frac{1}{r} \int_0^{+\infty} \sin(qr)\cdot \frac{q}{q^2+q_{TF}^2} \mathrm{d}q\\ =& \frac{-e}{4 \pi \varepsilon_0} \cdot \frac{1}{r} \cdot \frac{2}{\pi}\int_0^{+\infty} \sin(qr)\cdot \frac{q}{q^2+q_{TF}^2} \mathrm{d}q \end{align}

对

\begin{align} \frac{q}{q^2+q_{TF}^2} \end{align}

在 \(q \rightarrow +\infty\) 展开(与在 \(q_{TF}\rightarrow 0\) 时对 \(q_{TF}\) 展开相同):

import sympy as sym
q = sym.Symbol('q')
q_tf = sym.Symbol('q_TF')
s = sym.series(q/(q**2+q_tf**2),q,sym.oo,5)
print(sym.latex(s))

\frac{q_TF²}{q³} + \frac{1}{q} + O\left(\frac{1}{q⁵}; q→ ∞\right)

\begin{align} \frac{q}{q^2+q_{TF}^2} =\frac{q_{TF}^{2}}{q^{3}} + \frac{1}{q} + O\left(\frac{1}{q^{5}}; q\rightarrow \infty\right) \end{align}

计算 leading order

\begin{align} \int_0^{+\infty} \sin(qr)\cdot \frac{1}{q}\mathrm{d}q = \frac{1}{2i} \int_{-\infty}^{+\infty} e^{iqr}\cdot \frac{1}{q} \mathrm{d} q = \frac{1}{2i} \cdot \pi i \cdot 1 = \frac{\pi}{2} \end{align}

这正好是精确结果

\begin{align} \frac{-e}{4 \pi \varepsilon_0} \cdot \frac{e^{-q_{TF} r}}{r} \end{align}

在 \(q_{TF}\rightarrow 0\) 时对 \(q_{TF}\) 展开的结果的 leading order 相同.

问题1

Thomas-Fermi 近似结果与精确结果的 \(q_{TF}^2\) 的系数不同?

答: 因为展开后的结果中出现了发散项, 所以问题本身就不能展开.

问题2

为什么要对 \(q\) 在 \(q\rightarrow +\infty\) 展开? 积分的区间不是整个实轴吗?

答: 我之前对展开的目的理解有误. 展开的目的不是为了一阶一阶地计算积分. 展开目的是为了数值的计算积分. Thomas-Fermi 的积分本身可以精确解出, 没必要数值积分, 因此不必做展开. 展开的意义会在下面 RPA 的计算过程中体现.

RPA 的结果

RPA 给出

\begin{align} \varepsilon (\vec{q}) = 1+ \frac{q_{TF}^2}{q^2}g\left( \frac{q}{2k_F} \right) \end{align}

将上式代入 \(\phi (\vec{q})\) 有

\begin{equation} \varphi (\vec{q})= \frac{-e}{\varepsilon_0 V q^2 \left(1 +\frac{q_{TF}^2}{q^2}\right)} =\frac{-e}{\varepsilon_0 V \left(q^2 +q_{TF}^2 g\left( \frac{q}{2k_F} \right) \right)} \end{equation}

其中

\begin{align} g(u) = \frac{1}{2} \left( 1+\frac{1}{2u}(1-u^2)\ln \left| \frac{1+u}{1-u} \right| \right) \end{align}

对 \(\phi (\vec{q})\) 作 Fourier Transform

\begin{align} \varphi(\vec{r}) =\frac{-e}{2 \pi^2 \varepsilon_0} \cdot \frac{1}{r} \int_0^{+\infty} \sin(qr)\cdot \frac{q}{q^2+q_{TF}^2g( \frac{q}{2k_F} )} \mathrm{d}q \end{align}

级数展开

对

\begin{align} \frac{q}{q^2+q_{TF}^2g( \frac{q}{2k_F} )} \end{align}

在 \(q\rightarrow +\infty\) 展开

import sympy as sym
q = sym.Symbol('q')
q_tf = sym.Symbol('q_TF')
kf = sym.Symbol('k_F')
u = q/(2*kf)
*g=1
g = sym.Rational(1,2)*( 1+(1-u**2)/(2*u)*sym.log((u+1)/(u-1) ) )
s = sym.series(q/(q**2+q_tf**2)*g,q,sym.oo,10)
print(sym.latex(s))

\frac{\frac{256 k_F⁸}{63} - \frac{64 k_F⁶ q_TF²}{35} + \frac{16 k_F⁴ q_TF⁴}{15}}{q⁹} + \frac{\frac{64 k_F⁶}{35} - \frac{16 k_F⁴ q_TF²}{15} + \frac{4 k_F² q_TF⁴}{3}}{q⁷} + \frac{\frac{16 k_F⁴}{15} - \frac{4 k_F² q_TF²}{3}}{q⁵} + \frac{4 k_F²}{3 q³} + O\left(\frac{1}{q¹⁰}; q→ ∞\right)

\begin{align} \frac{q}{q^2+q_{TF}^2g( \frac{q}{2k_F} )} = \frac{\frac{256 k_{F}^{8}}{63} - \frac{64 k_{F}^{6} q_{TF}^{2}}{35} + \frac{16 k_{F}^{4} q_{TF}^{4}}{15}}{q^{9}} + \frac{\frac{64 k_{F}^{6}}{35} - \frac{16 k_{F}^{4} q_{TF}^{2}}{15} + \frac{4 k_{F}^{2} q_{TF}^{4}}{3}}{q^{7}} + \frac{\frac{16 k_{F}^{4}}{15} - \frac{4 k_{F}^{2} q_{TF}^{2}}{3}}{q^{5}} + \frac{4 k_{F}^{2}}{3 q^{3}} + O\left(\frac{1}{q^{10}}; q\rightarrow \infty\right) \end{align}

问题3: 接下来该怎么做? 所有的展开项代入积分都是发散的.

答: 展开算错了. 程序中把公式抄错了. 看来最好还是手算 check 一下比较好, 不然不会出现这么低级的错误. 级数展开的计算能力还有待提高…

正确的级数展开

import sympy as sym
q = sym.Symbol('q')
q_tf = sym.Symbol('q_TF')
kf = sym.Symbol('k_F')
u = q/(2*kf)
g = sym.Rational(1,2)*( 1+(1-u**2)/(2*u)*sym.log(sym.Abs((u+1)/(u-1)) ) )
s = sym.series(q/(q**2+q_tf**2*g),q,sym.oo,10)
print(sym.latex(s))

\frac{1}{q} - \frac{4 k_F² q_TF²}{3 q⁵} - \frac{16 k_F⁴ q_TF²}{15 q⁷} - \frac{64 k_F⁶ q_TF²}{35 q⁹} + O\left(\frac{1}{q¹⁰}; q→ ∞\right)

\begin{align} \frac{q}{q^2+q_{TF}^2g( \frac{q}{2k_F} )} = \frac{1}{q} - \frac{4 k_{F}^{2} q_{TF}^{2}}{3 q^{5}} - \frac{16 k_{F}^{4} q_{TF}^{2}}{15 q^{7}} - \frac{64 k_{F}^{6} q_{TF}^{2}}{35 q^{9}} + O\left(\frac{1}{q^{10}}; q\rightarrow \infty\right) \end{align}

数值计算 RPA 的结果

问题分析

首先, 积分

\begin{align} \varphi(\vec{r}) =\frac{-e}{2 \pi^2 \varepsilon_0} \cdot \frac{1}{r} \int_0^{+\infty} \sin(qr)\cdot \frac{q}{q^2+q_{TF}^2g( \frac{q}{2k_F} )} \mathrm{d}q \end{align}

难以解析地计算, 因此想要数值地计算.

其次, 被积函数收敛不够快.

经过前面对被积函数级数展开的分析, 发现被积分函数在远处的 leading order 为

\begin{equation} \label{eq:leading} \frac{\sin(qr)}{q} \end{equation}

也就是说, 被积函数在远处的贡献主要来自于 leading order 项 (\ref{eq:leading}), 而且这一项是可以解析地计算出结果的. 如果我们从被函数中把这项减去, 那么被积函数在远处就近似为 \(0\) 了, 也就是说积分会更快地收敛. 这就解决了积分收敛不够快的问题.

构造新的被积函数

构造函数

\begin{equation} F(q) = \sin(qr)\cdot \frac{q}{q^2+q_{TF}^2g( \frac{q}{2k_F} )} - \frac{\sin(qr)}{q} \end{equation}

新的 \(F(q)\) 函数会很快地收敛.

原来的积分就可以分解为

\begin{align} \varphi(\vec{r}) =&\frac{-e}{2 \pi^2 \varepsilon_0} \cdot \frac{1}{r} \int_0^{+\infty} \sin(qr)\cdot \frac{q}{q^2+q_{TF}^2g( \frac{q}{2k_F} )} \mathrm{d}q \\ =& \frac{-e}{2 \pi^2 \varepsilon_0} \cdot \frac{1}{r} \int_0^{+\infty} \left( F(q) + \frac{\sin(qr)}{q} \right)\mathrm{d}q \\ =& \frac{-e}{2 \pi^2 \varepsilon_0} \cdot \frac{1}{r} \left( \int_0^{+\infty} F(q) \mathrm{d}q + \frac{\pi}{2} \right) \end{align}

数值积分

下面数值地计算

\begin{align} \int_0^{+\infty} F(q) \mathrm{d}q \end{align}

python 程序如下

#导入数值计算包, 画图包, 积分包
import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate

#设置物理参量
qtf = 3/4
kf = 1/2

#数值积分的区间:(0,end)
end = 5

#定义函数 phi_r, g(u), phi_q, lo(q)是phi_q 的 leading order, F(q)是从 phi(q) 里边减去 leading order
def phi_r(r):
    def g(u):
	return 1/2*(1+1/(2*u)*(1-u**2)*np.log(np.abs((1+u)/(1-u))))
    def phi(q):
	phi = np.sin(q*r)*q/(q**2+qtf**2*g(q/(2*kf)))
	return phi
    def lo(q):
	lo = np.sin(q*r)/q
	return lo
    def F(q):
	F = phi(q) - lo(q)
	return F
    (Fres,err) = integrate.quad(F,0,end)
    phir = Fres + np.pi/2
    return phir


#计算N个不同的 r 的取值的数值积分 
N = 100


#画六个图中的第一个图
plt.subplot(321) #(321)的意思是整个图是3行2列,这个图画在第1个位置
a = 0.1 #画图的区间为[a,b],这里第一个图为 r 取 [0,10]
b = a+10
r = np.linspace(a,b,N)

#进行 N 次数值积分, 然后把 N 个点的数值积分的结果保存在数组 s 中
s = np.linspace(0,0,N)
for i in range(N):
    s[i-1] = phi_r(r[i-1])
#画图, 指定线的颜色和图例
plt.plot(r, s, color="red", label="Numerical Result")
#我们想要将数值的结果和 Friedel Oscillations 的结果做比较, 但它们之间会差一个常数k, 现在把这个常数 k 取 N 次数值积分中第一个点的比值
k = phi_r(a)/(1/(a**3)*np.cos(2*kf*a))
fo = k*1/(r**3)*np.cos(2*kf*r) #这是 Friedel Oscillations 的结果
plt.plot(r,fo, "--", color="blue", label="Friedel Oscillation") #画出 Friedel Oscillations 的结果

plt.xlabel("$r$") #设置横坐标的物理量
plt.ylabel("$\phi(r)$") #设置纵坐标的物理量

#以下三行是为了让坐标的数值以科学计数法显示
ax = plt.gca()  # 获取当前图像的坐标轴信息
ax.yaxis.get_major_formatter().set_powerlimits((0,1)) # 将坐标轴的base number设置为一位。
plt.legend()



#与画第一个图同样的方法, 画出剩下的5个图, 它们的区别在于 r 的区间不同
plt.subplot(322)
a = 10
b = a+10
r = np.linspace(a,b,N)

s = np.linspace(0,0,N)
for i in range(N):
    s[i-1] = phi_r(r[i-1])
plt.plot(r, s, color="red", label="Numerical Result")

k = phi_r(a)/(1/(a**3)*np.cos(2*kf*a))
fo = k*1/(r**3)*np.cos(2*kf*r)
plt.plot(r,fo, "--", color="blue", label="Friedel Oscillation")

plt.xlabel("$r$")
plt.ylabel("$\phi(r)$")

ax = plt.gca()  # 获取当前图像的坐标轴信息
ax.yaxis.get_major_formatter().set_powerlimits((0,1)) # 将坐标轴的base number设置为一位。

plt.subplot(323)
a = 20
b = a+10
r = np.linspace(a,b,N)

s = np.linspace(0,0,N)
for i in range(N):
    s[i-1] = phi_r(r[i-1])
plt.plot(r, s, color="red", label="Numerical Result")

k = phi_r(a)/(1/(a**3)*np.cos(2*kf*a))
fo = k*1/(r**3)*np.cos(2*kf*r)
plt.plot(r,fo, "--", color="blue", label="Friedel Oscillation")

plt.xlabel("$r$")
plt.ylabel("$\phi(r)$")

ax = plt.gca()  # 获取当前图像的坐标轴信息
ax.yaxis.get_major_formatter().set_powerlimits((0,1)) # 将坐标轴的base number设置为一位。


plt.subplot(324)
a = 50
b = a+10
r = np.linspace(a,b,N)

s = np.linspace(0,0,N)
for i in range(N):
    s[i-1] = phi_r(r[i-1])
plt.plot(r, s, color="red", label="Numerical Result")

k = phi_r(a)/(1/(a**3)*np.cos(2*kf*a))
fo = k*1/(r**3)*np.cos(2*kf*r)
plt.plot(r,fo, "--", color="blue", label="Friedel Oscillation")

plt.xlabel("$r$")
plt.ylabel("$\phi(r)$")

ax = plt.gca()  # 获取当前图像的坐标轴信息
ax.yaxis.get_major_formatter().set_powerlimits((0,1)) # 将坐标轴的base number设置为一位。


plt.subplot(325)
a = 100
b = a+10
r = np.linspace(a,b,N)

s = np.linspace(0,0,N)
for i in range(N):
    s[i-1] = phi_r(r[i-1])
plt.plot(r, s, color="red", label="Numerical Result")

k = phi_r(a)/(1/(a**3)*np.cos(2*kf*a))
fo = k*1/(r**3)*np.cos(2*kf*r)
plt.plot(r,fo, "--", color="blue", label="Friedel Oscillation")

plt.xlabel("$r$")
plt.ylabel("$\phi(r)$")

ax = plt.gca()  # 获取当前图像的坐标轴信息
ax.yaxis.get_major_formatter().set_powerlimits((0,1)) # 将坐标轴的base number设置为一位。


plt.subplot(326)
a = 200
b = a+10
r = np.linspace(a,b,N)

s = np.linspace(0,0,N)
for i in range(N):
    s[i-1] = phi_r(r[i-1])
plt.plot(r, s, color="red", label="Numerical Result")

k = phi_r(a)/(1/(a**3)*np.cos(2*kf*a))
fo = k*1/(r**3)*np.cos(2*kf*r)
plt.plot(r,fo, "--", color="blue", label="Friedel Oscillation")

plt.xlabel("$r$")
plt.ylabel("$\phi(r)$")

ax = plt.gca()  # 获取当前图像的坐标轴信息
ax.yaxis.get_major_formatter().set_powerlimits((0,1)) # 将坐标轴的base number设置为一位。



plt.suptitle("The Results of RPA") #画出整个图的标题
plt.show() #显示图片

结果对比

将数值结果与 Friedel Oscillations 的结果

\begin{align} \varphi(\vec{r}) \sim \frac{1}{r^3} \cos (2 k_fr) \end{align}

进行对比

可以看出, Friedel Oscillations 的结果在 \(r\) 较大时的结果比较好.

参考文献

Friedel Oscillation 的原始文献 the shielding of a fixed charge in a high-density electron gas http://www.doc88.com/p-9512851691956.html

Wolfgang Nolting, Fundamentals of Many-body Physics

致谢

感谢导师 Ran Qi 的指导

感谢 Fan Yang 师兄的讨论